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Abstract 

In  1993,  Scherpen  generalized  the  balanced  truncation 
method  to  the  nonlinear  setting.  However,  the  Scherpen  pro¬ 
cedure  is  not  easily  computable  and  has  not  yet  been  applied 
in  practice.  We  offer  methods,  tools,  and  algorithms  for  com¬ 
puting  the  energy  functions  and  coordinate  transformations 
involved  in  the  Scherpen  theory  and  procedure  for  nonlin¬ 
ear  balancing.  We  then  apply  our  approach  to  derive,  for  the 
first  time,  balanced  representations  of  nonlinear  state-space 
models. 

1  Introduction 


Definition  1.1  (Controllability  Function)  The  controlla¬ 
bility  function,  Lc  :  R”  — >  R ,  M  system  (l)-(2)  is  defined 
by 


Lc  {x0)  — 


min 

u  G  L2(—oo,  0) 
x(—oo)  =  0  ,  x(0)  =  Xq 


IN*)  II2  dt 


(3) 

□ 


Definition  1.2  (Observability  Function)  The  observability 
function,  La  :  Rn  -y  R ,  for  system  (l)-(2)  is  defined  by 


L0  (To)  —  2 
®(0)  =  Xq  , 


IN*)  II2  dt 


10 

u(t)  =  0,  t  >  0 


(4) 


□ 


This  paper  addresses  the  problem  of  computability  pertain¬ 
ing  to  the  Scherpen  theory  and  procedure  for  balancing  of 
nonlinear  systems  [1,  2].  We  offer  methods,  tools,  and  algo¬ 
rithms  toward  computing  balanced  realizations  for  asymptot¬ 
ically  stable  affine  nonlinear  control  systems,  i.e.,  state-space 
models  of  the  form 


We  will  use  the  following  results  pertaining  to  the  control¬ 
lability  function  of  a  nonlinear  system. 

Proposition  1.3  In  the  case  of  an  linear  time-invariant 
(LTI)  system  the  controllability  function  specializes  to  the 
quadratic  form 


x{t)  =  f  {x(t))  +  ^2  gi(x(t))  lift)  (1) 

i=  1 

y(t)  =  h(x(t))  (2) 

where  u  =  (u  i, . . . ,  um)  G  U  C  R”\  y  =  (yi,.--,yp)  G 
Rp,  and  x  =  (aq , . . . ,  xn)  are  local  coordinates  for  a  smooth 
state-space  manifold  M.  The  maps  f,g i , . . . ,  gm  are  smooth 
and  we  assume  that  /( 0)  =  0  and  h( 0)  =  0. 

We  say  that  /  is  stable  (asymptotically  stable)  if  0  is  a 
stable  (asymptotically  stable)  equilibrium  for  x  =  /(.:/:),  and 
normally  assume  asymptotic  stability  of  /.  We  refer  to  the 
triple  (/,  g,  h )  as  a  realization  of  the  nonlinear  system. 

The  two  main  objects  that  we  work  with  in  nonlinear  bal¬ 
ancing  are  the  controllability  and  observability  energy  func¬ 
tions  of  the  system,  defined  as  follows. 


Lc(x0)  =  ^x0TWc  1  Xq  (5) 

where  the  symmetric  positive-definite  matrix  Wc  is  the  famil¬ 
iar  controllability  Gramian.  □ 


Theorem  1.4  (Scherpen  [1,  2])  Suppose  that  0  is  an  asymp¬ 
totically  stable  equilibrium  of  —  +  g  gJ  ^  on  a 

neighborhood  V  of  0.  Then,  for  all  x  G  V,  Lc  is  the  unique 
smooth  solution  of 


°  =  ~^x  ^  ^  +  \  ~^x  ^  9  ^  ^  ^ 


(6) 


with  supplementary  condition  Lc(0)  =  0  and  under  the  as¬ 
sumption  that  (6)  has  a  smooth  solution  on  V. 


In  this  paper  we  will  consider  two  of  the  fundamental 
problems  involved  in  computing  balanced  realizations  for 
nonlinear  systems.  First,  we  consider  the  problem  of  com¬ 
puting  the  controllability  energy  function  without  solving 
the  family  of  optimal  control  problems  implied  in  its  def¬ 
inition,  or  solving  the  associated  Hamilton-Jacobi-Bellman 
(HJB)  equation  (6).  Then,  we  consider  the  problem  of  com¬ 
puting  a  Morse  coordinate  transformation  under  which  the 
controllability  function  is  locally  quadratic  on  a  neighbor¬ 
hood  of  0.  We  apply  the  methods  developed  in  this  paper  to 
two  example  systems:  a  forced  damped  pendulum  system, 
and  a  forced  damped  double  pendulum  system. 

Methods  for  computing  the  observability  function  and  the 
balancing  coordinate  transformation,  algorithmic  details,  and 
additional  numerical  studies  can  be  found  in  [3].  We  note 
here  that  the  authors  have  developed  a  MATLAB  toolbox  for 
implementing  the  methods  described  herein,  which  was  used 
for  performing  all  of  the  numerical  studies. 

2  Stochastic  Methods  for  Computation 

We  seek  a  method  for  computing  the  controllability  energy 
function  without  solving  the  family  of  optimal  control  prob¬ 
lems  implied  in  its  definition,  or  solving  the  associated  HIB 
equation  (6).  In  this  section  we  present  an  approach,  based 
primarily  on  the  theory  of  stochastically  excited  dynamical 
systems  (see,  e.g.,  [4,  5]),  for  computing  an  estimate  of  the 
controllability  function.  We  show  that  in  certain  situations 
the  method  provides  an  exact  solution.  Below,  for  integer  k, 
we  let  k  denote  the  set  {1, . . . ,  k}. 

By  a  stochastically  excited  dynamical  system  we  mean  a 
control  system  for  which  the  m  components  of  the  input,  u,, 
i  £  m,  have  been  replaced  by  the  sample  paths  of  m  Gaus¬ 
sian  white  noises,  {((t)i,  t  G  R+},  i  G  m.  The  state  equa¬ 
tion  is  given  by 

,  m 

-Xt=f(Xt)  +  '£gi(Xt)(Ct)i  (7) 

i—  1 

The  white  noise  driven  system  (7)  is  interpreted  correctly 
via  the  stochastic  differential  equation  (SDE)  (given  elemen¬ 
twise)  for  i  £  n 


{dXt)i  =  fi  {Xt)  dt  +  Y^gi  (Xt)  (■ dWt ) .  (8) 

i= 1 

where  for  i  E  n 


1  n  m  rx 

fi  W+2EE9I  (**)  9jk  (Xt) 

j—1  k—  1  ^ 

(9) 

includes  the  so  called  correction  term,  gij  =  (g:j) .,  i  £ 
n,  j  £  m,  and  where  solutions  to  (8)  are  defined  in  terms 
of  a  corresponding  stochastic  integral. 

The  state  Xt  is  a  Markov  process  with  transition  proba¬ 
bility  density  p(x,  t;  y,  s).  Time  evolution  of  p(x,  t;  y,  s )  is 


fi  (Xt)  = 


governed  by  the  Fokker-Planck  equation,  given  by 
dp 


dt 


(x,t;y,s)  = 
n  d 


i=  1 


-y  n  n  q2 

(hj  (x)p(x,i;y,s))  (10) 


2  -4-'  ^  dxi  dx 

»= 1  j=  1 

with  initial  condition 

P(x,s;y,s)  =  S(x-y) 


and  where 


bij  ( x )  =  9ik  gik  ^  = 

k=  1 

Finally,  in  the  steady-state.  Equation  (10)  simplifies  to  the 
stationary  Fokker-Planck  equation 


[s0*0]  bWlT 


0 


_E  (f*  (x)  p°°  (x))  (H) 

2—1 

n  n  r,2 

+2  EEss;i*»  (*)!>»(*)) 

*= 1  j= 1  J 


where  (x)  denotes  the  stationary  probability  density  (if  it 
exists). 

In  the  LTI  case,  the  transition  density  function  p(x,  t;  y,  s) 
describing  the  random  properties  of  the  state  process  Xt  is 
Gaussian,  and  the  stationary  density  p ^  has  zero  mean  and 
covariance  equal  to  the  controllability  Gramian  matrix  Wc 
(see,  e.g.,  [6]),  i.e., 

Poo(x)  =  [(27r)n  det(Wc)]-1/2  exp  ^~^xJ  W^1  x^j 

(12) 

Thus,  in  the  LTI  case,  the  controllability  function  Lc  and 
the  stationary  density  p^  are  related  exactly  by 

Poo(x)  =  [(27 r)"  det(Wc)p1/2  exp(-Lc(x))  (13) 


and 

Lc(x)  =  -log  (poc(x))  +  log  ([(2tt)”  det  (Wc)]~1,2^j 

(14) 

In  the  nonlinear  setting,  the  density  p[x,  t:  y,  s),  and  in 
particular  the  stationary  density  Poo(x),  are  not,  in  general, 
Gaussian,  nor  determined  completely  by  their  mean  and  co- 
variance,  i.e.,  higher  order  moments  are  involved.  However, 
because  the  balancing  coordinate  transformation  is  local  to  a 
neighborhood  of  the  origin,  we  are  mainly  interested  in  cap¬ 
turing  a  local  characterization  of  the  controllability  function. 
In  light  of  this.  Equation  (14)  suggests  that  a  useful  approxi¬ 
mation  of  Lc  is  defined  by 

L'c(x)  =  -log  (Poo(x))  +  C 


(15) 


where  C  is  a  normalizing  constant,  dependent  on  the  partic¬ 
ular  system,  such  that  L'c  (0)  =  0.  By  Equation  (14),  L'c 
specializes  to  the  exact  Lc  in  the  LTI  case. 

There  exist  certain  nonlinear  systems  for  which  Equa¬ 
tion  (15)  provides  an  exact,  rather  than  approximate,  formula 
for  the  controllability  function.  As  one  simple  example,  con¬ 
sider  the  Langevin  process  Xt  governed  by  the  first-order 
SDE 

dXt  = -X  (j)(Xt)  +  dWt  (16) 

where  <f>  :  R"  — >■  R"  is  a  Cl  map  such  that  —  V  (f>  is  asymp¬ 
totically  stable.  The  Maxwell-Boltzmann  density 

PooB(a;)  =  C  exp  (—2  <t>(x))  (17) 


satisfies  the  stationary  Fokker-Planck  equation  where  C  is  a 
constant  such  that  f  p^B  =  1.  Now,  using  (15),  define 

LfB  (*)  =  -log  (p™B(z))  +  log  (C)  =  2  (18) 


which  is  the  unique  controllability  energy  function  for  the 
Langevin  system,  i.e.,  stable  affine  nonlinear  system  with 
f(x)  =  —V  <j>(x)  and  g(x)  =  I. 

Now  we  seek  conditions  under  which  a  broader  class  of 
systems  admits  an  exact  relationship  between  the  stationary 
density  and  the  controllability  function.  In  particular,  we 
consider  second-order  mechanical  systems  with  a  Hamilto¬ 
nian  structure  perturbed  by  dissipation  and  forcing.  We  de¬ 
termine  conditions  under  which  the  controllability  function 
for  such  a  system  can  be  expressed  exactly  in  terms  of  the  sta¬ 
tionary  density  for  the  corresponding  stochastically  excited 
system.  We  adopt  and  modify  somewhat  the  notation  and 
framework  of  Fuller  [7]  and  Zhu  and  Yang  [8],  These  au¬ 
thors  have  presented  conditions  under  which  exact  solutions 
of  the  stationary  Fokker-Planck  equation  can  be  derived.  We 
show  that  in  certain  cases,  the  same  conditions  are  sufficient 
for  expressing  the  controllability  function  in  terms  of  the  sta¬ 
tionary  density,  while  in  other  cases,  additional  conditions 
are  required. 

We  consider  a  forced,  dissipatively  perturbed,  n-DOF 
Hamiltonian  system.  Let  q  =  (qi,  ■  ■  ■  ,qn)  £  R”  and 
p  =  (pi, ,  Pn)  £  R"  denote,  respectively,  the  gener¬ 
alized  displacements  and  momenta.  Let  the  Hamiltonian 
H'  =  H'  ( q,p ),  i.e.,  the  sum  of  the  kinetic  and  potential  en¬ 
ergies  of  the  system,  be  C 2.  Let  c';  =  cF  (q,p)  for  i,j  £  n 
be  C 1  functions  representing  nonlinear  dissipation  coeffi¬ 
cients.  Let  dij  =  dij  ( q,p )  for  i,j  £  n  be  C2.  The  system 
that  we  consider  is  governed  by  the  equations  of  motion,  for 
i  £  n 


Qi 


Pi 


dH' 

dpi 

dH' 

dqi 


dH' 

dPj 


m 

+ dik  u 

k= i 


(19) 

(20) 


The  system  is  realized  in  standard  state-space  form  with 
coordinates  x  =  ( q,p )  £  R2™  and 


dH' 

dpi 


fi 


C 9k)  i 
C 9k)  i 


dH' 

dqi 


n 


E  cij 


dH' 

dPj 


i  =  n  +  1, . . . ,  2n 


0  k=l,...,m  (22) 

dik  i  =  n  +  1, . . . , 2n;  k  =  l,...,m 


The  output  map  h  is  irrelevant  for  purposes  of  the  discussion 
here. 

The  corresponding  stochastically  excited  system  is  gov¬ 
erned  by  the  SDEs,  for  i  £  n 


dQi 


dH' 

m 


dt 


(23) 


dPi 


( OH'  ^  ,  dH_ 
[  <)Q;  2^  r"  nr, 

\  3=1 


ddik 

dPj 


(24) 

dt 


+  ^2dik  (dWt). 

k=  1 


where  we  have  adopted  the  usual  notation  by  substituting 
Q  for  q  and  P  for  p  when  dealing  with  the  corresponding 
random  variables. 

It  is  usually  the  case  that  the  correction  terms 


1  n  m  r\  i 
1  v  OCLik 


djk ,  i  £  n,  can  be  split  into  two  parts:  one 


2  op 

j= l  fc= l  J 

which  modifies  the  conservative  forces  and  the  other  that 


modifies  the  damping  forces  (see  [8]).  This  allows  Equa¬ 
tions  (23)  and  (24)  to  be  rewritten,  for  i  £  n 


dQi 


dPi 


+  Y,dik  (dWt\ 

fc= l 


dt 


(25) 

(26) 


for  appropriately  defined  H  and  ct] . 

The  stationary  Fokker-Planck  equation  governing  the  sta¬ 
tionary  transition  density  p^  =  px  ( q,p )  associated  with 
the  SDEs  (25)  and  (26)  is  given  by 


(27) 


-  E 

i= 1 
n 

+  E 


d 

(dH  n  ^ 

.  9  I 

( dH  W 

dQi 

\dPiPo°) 

]  +  m  1 

laQ-pooJJ 

o 

dPi 


E 

VJ=1 


dH_ 


cn  dp  P. 


n 


d2 


3  =  1 


dPi  dP i 


(pij  Po 


where 


bij  =  ^2  djk  =  [d]  [d]1 


k= 1 


and  subject  to  boundary  conditions  (vanishing  probability 
flow) 

,  dH 

hm  — —  Poo  =  0 

II  (g>p)  dPi 


fi 


i  =  1, . . . ,  n 


(21) 


(28) 


and 


Observe  that  the  first  summation  term  on  the  right-hand- 
side  of  (27)  is  equal  to  the  Poisson  bracket  of  p-^  and  H, 
i.e., 

{Poc,H}  (30) 

dH  dpoo  dH_  dpao 

dPi  dQi  dQi  dPi 

' _ d_  (dH  \  d_  (dH_  V 

.dQi  \dPiPoo)  +  dPi  \dQPo°). 

Thus,  we  can  rewrite  the  stationary  Fokker-Planck  equa¬ 
tion  (27)  as 


=  E 

i—  1 
n 

=  E 


Remark  2.4  The  equipartition  of  energy  condition  imposes 
a  severe  restriction  on  the  class  of  systems  for  which  (33)  is 
the  stationary  density.  □ 


The  controllability  function  Lc  uniquely  satisfies  the  HJB 
equation  (6),  which,  for  realization  (/,  g)  given  in  Equa¬ 
tions  (21)  and  (22),  takes  the  form 


0  =  {Lc,H}  +  Y, 

i= 1 


dLc  ”  / 

'  dHl 

0LC\ 

F§< 

v  Cijdfj+2  ” 

dpj 

(34) 

The  relationship  between  the  stationary  density  and  the  con¬ 
trollability  function,  and  an  exact  formula  for  the  latter  in 
terms  of  the  Hamiltonian,  are  given  in  the  following  result. 


Theorem  2.5  Consider  the  forced,  dissipatively  perturbed, 
n-DOF  Hamiltonian  control  system,  governed  by  the  evolu¬ 
tion  equations  (19)  and  (20),  and  realized  by  ( f,g )  given  in 
Equations  (21 )  and  (22).  Under  the  conditions  stated  in  The¬ 
orem  2.1,  the  unique  controllability  energy  function  for  the 
system  is  given  by 


0 

=  {Poo,H} 


(3D  Lc  (q,p)  =  -log  (p^q,  p)) +  C' 

=  2  eH(q,p)  +  C'  (35) 


+ 


E 


3= 1 


d2 

OP i  dPj 


( pij  Poo  ) 


The  relationship  between  the  stationary  density  and  the  con¬ 
trollability  function,  and  an  exact  formula  for  the  latter  in 
terms  of  the  Hamiltonian,  are  given  in  the  following  results. 
We  first  consider  the  case  where  the  parameters  cC  and  dik 
are  independent  of  q  and  p,  i.e.,  are  constants.  In  this  situa¬ 
tion,  the  correction  term  vanishes,  H  =  H',  and  Cy  =  c'? . 
The  following  is  modified  from  Fuller  [7]. 


Theorem  2.1  (Fuller  [7])  Consider  the  stochastically  ex¬ 
cited  system  corresponding  to  the  forced,  dissipatively 
perturbed,  n-DOF  Hamiltonian  system  governed  by  the 
SDEs  (25)-(26)  where  H  is  the  Hamiltonian.  Suppose  that 
the  coefficients  Cij  and  dik  are  independent  of  q  and  p.  Fur¬ 
thermore,  suppose  that  the  following  constant  ratio  holds  for 
all  i,j  £  n: 

-EL  =  l  =  constant  (32) 

bij 

Then  the  unique  stationary  density  p^  that  satisfies  Equa¬ 
tion  (31)  is 


where  is  the  stationary  density  of  the  corresponding 
stochastically  excited  system  and  C'  is  a  constant  such  that 
Lc(  0,0)  =  6. 

Proof  It  is  necessary  and  sufficient  to  show  that  Lc  satisfies 
Equation  (34).  Since  Lc  is  a  functional  of  H ,  we  have  that 

dL  dH 

{LC,H}  =  0.  Furthermore,  —L  =  21  — — ,  so  that  Equa- 

dPj  ac¬ 

tion  (34)  becomes 


0  = 


E 

i= 1 
n 

E 


dH^dH 

dPi  ^  dpj  13 


(37) 


which  is  clearly  satisfied  given  the  equipartition  of  energy 
condition  (32).  ■ 

We  now  consider  the  more  general  situation  where  the  pa¬ 
rameters  c';  and  dik  are  permitted  to  be  functions  of  q  and p. 
The  following  is  modified  from  Zhu  and  Yang  [8]. 


Poo  (q,p)  =  Cexp(-2£H  (q,p))  (33) 

where  C  is  a  constant  such  that  f  p^  =  1. 

Remark  2.2  The  density  ( 33 )  is  of  Maxwell-Boltzmann 
form.  □ 


Theorem  2.6  (Zhu  and  Yang  [8])  Consider  the  stochasti¬ 
cally  excited  system  corresponding  to  the  forced,  dissipa¬ 
tively  perturbed,  n-DOF  Hamiltonian  system  governed  by 
the  SDEs  (25)-(26)  where  H  is  the  Hamiltonian.  Suppose 
that  the  following  ratio  holds  for  all  i  £  n  and  for  some 
functional  h  of  H: 


Remark  2.3  The  condition  (32)  is  referred  to  as  the  equipar¬ 
tition  of  energy  condition.  The  terminology  derives  from  the 
situation  in  statistical  mechanics  where  each  DOF  of  a  multi- 
particle  system  is  associated  with  the  same  mean  energy.  □ 
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=  h(H) 


(38) 


Then  the  unique  stationary  density  that  satisfies  Equa-  3  Computing  the  Morse  Coordinate  Transfer- 
tion  (31)  is  mation 


Poo  (q,p)  =  C  exp  J 


H(q,p) 


h(u)  duj  (39) 
where  C  is  a  constant  such  that  f  p^  =  1. 


Remark  2.7  The  density  (39)  is  of  Maxwell-Boltzmann 
form.  □ 


Remark  2.8  The  condition  (38)  is  analogous  to  an  equipar- 
tition  of  energy  condition,  again  imposing  a  severe  restric¬ 
tion  on  the  class  of  systems  for  which  (39)  is  the  stationary 
density.  □ 


The  relationship  between  the  stationary  density  and  the 
controllability  function,  and  an  exact  formula  for  the  latter  in 
terms  of  the  Hamiltonian,  are  given  in  the  following  result. 


Theorem  2.9  Consider  the  forced,  dissipatively  perturbed, 
n-DOF  Hamiltonian  control  system,  governed  by  the  evolu¬ 
tion  equations  (19)  and  (20),  and  realized  by  ( f,g )  given  in 
Equations  (21)  and  (22).  Suppose  that  the  following  ratio 
holds  for  all  i  £  n  and  for  some  functional  r  of  H: 
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dH 


dH 

dpj 


=  r(H) 


(40) 


Then  the  unique  controllability  energy  function  for  the  system 
is  given  by 


rH(q>p) 

Lc(q,p)  =  2  /  r(u)du  +  C' 

Jo 


(41) 


where  C'  is  a  constant  such  that  Lc  (0,  0)  =  0.  Furthermore, 
if  the  bij  are  independent  of  p  then 


Lc  ( q,P )  =  -log  (Poo(g,p))  +  C'  (42) 


where  p^  is  the  stationary  density  of  the  corresponding 
stochastically  excited  system. 

Proof  Assume  that  Lc  ( q,p )  =  <p  ( H(q,p ))  for  some  func¬ 
tional  (f>  of  H.  Then  {LC,H}  =  0.  The  HJB  equation  (34) 
can  be  written 


Recall  that  for  an  LTI  system,  the  energy  functions  Lc  and 
L0  globally  take  the  form  of  quadratic  functions  given,  for 
example,  by  (5).  We  wish  to  generalize  the  linear  balanc¬ 
ing  procedure  to  the  nonlinear  setting,  but  the  functions  Lc 
and  La  are  not,  in  general,  quadratic.  However,  we  can 
appeal  to  some  important  results  from  critical  point  theory 
(see,  e.g.,  [9])  in  order  to  find  a  change  of  coordinates  un¬ 
der  which  a  smooth  function  takes  a  quadratic  form  locally 
around  a  non-degenerate  critical  point.  The  key  result  is  the 
Morse  lemma  [10],  which  guarantees  the  existence  of  the 
desired  canonical  form  for  functions  with  a  non-degenerate 
critical  point  (called  Morse  functions)  defined  on  a  finite¬ 
dimensional  manifold,  and  an  analogous  result  ofPalais  [11], 
which  generalizes  the  notion  to  functions  defined  on  a  Hilbert 
space.  Our  presentation  of  the  established  results  are  based 
mainly  on  those  found  in  [9,  12], 

Theorem  3.1  (Morse-Palais)  Let  f  be  a  smooth  real-valued 
function  defined  on  an  open  neighborhood  O  of  0  in  the 
Hilbert  space  £.  Assume  that  /( 0)  =  0  and  that  0  is  a  non¬ 
degenerate  critical  point  of  f.  Then  there  exists  a  neighbor¬ 
hood  U  C  O  of  0,  a  local  change  of  coordinates  (j>  on  U,  and 
an  invertible  symmetric  operator  A  such  that 

f  (x)  =  (A4>(x),<i>(x))£  x£U  (44) 

□ 

Corollary  3.2  Let  f  be  a  smooth  real-valued  function  de¬ 
fined  on  an  open  neighborhood  O  of  0  in  the  Hilbert  space 
£.  Assume  that  /( 0)  =  0  and  that  0  is  a  non-degenerate 
critical  point  of  f.  Then  there  exists  a  neighborhood  U  C  O 
of  0,  a  local  change  of  coordinates  z  =  t;(x)  on  U,  and  an 
orthogonal  decomposition  £  =  T  -\-  T _L  such  that  if  we  write 
z  =  £ (x)  =  u  +  v  with  u  £  T  and  v  £  J-1-  then 

f(z)  =  f(£(x))  =  (u,u)£-{v,v)£  x£U  (45) 

□ 

Remark  3.3  Consider  the  special  case  where  £  =  R”  and 
critical  point  0  has  index  r.  Define  z  =  ^(x)  for  x  £  U  and 
ip  =  £-1  on  £((/).  Then  Corollary  3.2  implies  that 


dH  dfi 
dpi  dH 


(  dH_d±_0  dH\ 

[  dp,  dH  °ij  dp J 


(43) 


f(x)  =  f  W(z))  =  E  zi  (46> 

i—1 


dd> 

which  is  clearly  satisfied  if  we  assign  — —  =2  r  ( H )  where 

dH 

r  (H)  is  defined  by  Equation  (40).  Thus,  the  desired 
functional  (j)  is  obtained  through  integration  yielding  Equa¬ 
tion  (41).  Moreover,  if  the  bij  are  independent  of  p  then 
2  r  ( H )  =  h  ( H )  where  h  (H)  is  defined  in  Equation  (38). 
In  that  case  Equation  (42)  holds.  ■ 

Additional  details  appear  in  [3], 


In  the  new  coordinates,  the  function  f  is  said  to  be  in  spher¬ 
ical  quadratic  form.  The  transformation  is  illustrated  in  Fig¬ 
ure  1.  □ 


Definition  3.4  (Morse  Coordinate  Transformation)  A 

change  of  coordinates  ip  satisfying  (46)  is  said  to  be  a  Morse 
coordinate  transformation  for  f  around  0.  □ 


with  symmetry  property 


Figure  1:  An  example  of  a  Morse  function  on  R2  (with 
level  contours)  before  and  after  transformation  to  spherical 
quadratic  form. 


The  original  proof  by  Morse  uses  the  Gram-Schmidt  or- 
thogonalization  process  which  is  essentially  a  coordinate-by- 
coordinate  induction  argument.  The  generalization  by  Palais 
is  proved  without  a  coordinate-wise  procedure,  which  we  use 
to  our  advantage  for  purposes  of  computation.  For  purposes 
of  brevity,  the  proofs  do  not  appear  here.  We  proceed  directly 
to  the  algorithms. 

We  present  here,  somewhat  loosely,  an  algorithm  for  nu¬ 
merical  implementation  of  Theorem  3.1  and  Corollary  3.2. 
The  algorithm  is  presented  more  rigorously  in  [3]  where  a 
computational  framework  is  introduced.  The  algorithm  takes 
a  Morse  function  /  and  returns  a  neighborhood  U,  Morse  co¬ 
ordinate  transformation  <j>,  and  invertible  symmetric  matrix  A 
under  which  /  takes  the  desired  form  (44)  on  U.  An  addi¬ 
tional  algorithm  takes  d,  A,  and  U  and  returns  a  coordinate 
transformation  £  under  which  /  takes  the  spherical  quadratic 
form  (45).  The  main  building  blocks  are  as  follows. 


smooth  function  decomposition  Given  smooth  real-valued 
function  /  with  critical  point  at  0,  return  smooth  func¬ 
tions  <7i,  i  £  n,  such  that 


n 

fix)  =  y  g%{x)  Xi  i  6  b,  x  £  O  (47) 

i—1 


and 


i  £n 


(48) 


df 

1.  Compute  (approximate)  partial  derivatives  — — , 

OXi 


i  £  n. 

2.  For  each  point  x  in  the  domain  of  defi¬ 
nition  of  /,  compute  (approximate)  integrals 

9i{x)  = 


[  (tx)  dt ,  i  £  n. 
Jo  dxi 


Morse  function  decomposition  Given  Morse  function  /, 
i.e.,  /  has  non-degenerate  critical  point  at  0,  return 
smooth  functions  hij,  i,j  £  n  such  that 


n  n 

f(X)  =  YYl  hijix)xixJ 

i=l  3= 1 


i,j£n,  x  £  O 
(49) 


hij(x)  =  hji(x)  i,j  £  n,  x  £  O  (50) 


and  such  that 


(0) 


1  d2f 

2  dxi  dxj 


(0)  =  1d2/(0) 


(51) 


This  is  accomplished  via  n  + 1  smooth  function  decom¬ 
positions. 


1.  Apply  the  smooth  function  decomposition  to  / 
yielding  g^i  £n. 

2.  Apply  the  smooth  function  decomposition  to  each 
of  the  gi  yielding  hij,  i,j  £  n. 

matrix  square  root  Given  matrix  B  close  to  the  identity,  re¬ 
turn  its  square  root  C,  i.e.,  B  =  C2.  The  matrix  B  must 
satisfy 

!  H  -  B  ||  <  1  (52) 

In  that  case,  the  following  algorithm  converges  to  a 
fixed  point  corresponding  to  the  desired  matrix  C  = 

B1'2. 

Ck+i  =  Ck  +  \{B~C2)  k  =  0,1,.. (.53) 

Co  =  I 


The  convergence  of  the  sequence  {Ck  \  to  the  fixed 
point  f?1/2  can  be  shown  to  be  a  consequence  of  the 
contraction  mapping  principle. 

Morse-Palais  transformation  Given  Morse  function  /,  re¬ 
turn  neighborhood  U,  coordinate  transformation  </>,  and 
invertible  symmetric  matrix  A  such  that  Equation  (44) 
holds. 


1.  Apply  the  Morse  function  decomposition  to  / 
yielding  hij,  i,j  £  n.  Let  H  ( x )  =  [hij  (x)]  and 
A  =  H  ( 0).  Note  that  the  nondegeneracy  of  0  en¬ 
sures  the  invertibility  of  A. 

2.  For  each  point  x  in  the  domain  of  definition  of  /: 

(a)  Compute  the  solution  B  of  the  matrix  equa¬ 
tion  AB  =  H  ( x ).  Note  that  B  =  II  at  x  =  0. 

(b)  If  ||  H  —  i?  ||  <  1  then: 

i.  Apply  the  matrix  square  root  algorithm  to 
compute  C  =  B1/2. 

ii.  Let  (f)  (x)  =  C  x. 

iii.  Include  the  point  x  in  the  neighborhood 
U. 

(c)  Otherwise,  do  not  include  the  point  x  in  the 
neighborhood  U  and  no  further  calculations 
apply. 

This  procedure  provides  an  estimate  of  the  neighbor¬ 
hood  U  for  which  the  function  can  be  transformed  to  the 
canonical  quadratic  form.  It  is  possible  that  the  maximal 
neighborhood  is  larger. 


spherical  transformation  Given  transformation  </>  and  in¬ 
vertible  symmetric  matrix  A  such  that  Equation  (44) 
holds,  return  index  r  and  coordinate  transformation  ip 
such  that  Equation  (46)  holds. 

1 .  Compute  the  spectral  decomposition  of  matrix  A, 
i.e.,  A  =  V  AVt. 

2.  Let  E  =  diag(|Ai| , . . . ,  |A„|). 

3.  Let  r  equal  the  number  of  A,  such  that  A,;  <  0. 

4.  Let  R  =  EVJ. 

5.  Lor  each  point  x  in  the  domain  of  definition  of  /, 

let  ip  (x)  =  Rcp{x). 

Remark  3.5  The  terminology  “Morse  function  decomposi¬ 
tion"  is  somewhat  misleading  since  the  decomposition  (49) 
merely  requires  a  critical  point  that  is  not  necessarily  non¬ 
degenerate.  However,  we  adopt  the  terminology  for  lack  of  a 
better  name  and  because  we  are  applying  the  decomposition 
to  Morse  functions.  □ 


Ligure  2:  Planar  pendulum  system  with  massless  shaft,  lin¬ 
ear  torsional  damping,  linear  torsional  stiffness,  and  torque 
input  applied  at  the  rotary  joint.  Parameter  values  used  for 
numerical  studies:  m  =  1/40,  b  =  2,  k  =  1,  L  =  20, 
G  =  10. 


4  Applications 

In  this  section  we  illustrate  our  methods  and  algorithms  by 
applying  them  to  two  examples  of  rigid  link  mechanical 
systems.  We  compute  a  balanced  realization  for  a  forced 
damped  pendulum  system,  and  take  steps  toward  balancing  a 
forced  damped  double  pendulum  system. 

The  first  example  that  we  consider  is  a  simple  pendu¬ 
lum  system  as  illustrated  in  Ligure  2.  The  system  incorpo¬ 
rates  linear  torsional  damping,  linear  torsional  stiffness,  and 
a  torque  input  at  the  rotary  joint.  We  assume  that  the  shaft  is 
massless  and  that  the  pendulum  moves  only  in  the  plane.  We 
consider  the  case  where  the  joint  angle  is  measured  (position 
read-out).  (See  [3]  for  the  case  of  velocity  read-out.) 

It  is  beneficial  to  study  the  pendulum  as  an  example  be¬ 
cause 

•  it  is  nearly  linear,  so  we  can  use  linear  theory  to  obtain  a 
good  estimate  of  the  correct  results  for  comparison;  and 

•  in  previous  sections  we  have  studied  second-order  me¬ 
chanical  systems  and  obtained  an  exact  formula  for  the 
controllability  function. 

We  obtain  a  state-space  realization  (/,  <j.  h)  for  the  pendu¬ 
lum  system  via  Lagrangian  mechanics.  Let  the  generalized 
position  q  and  velocity  q  correspond  to  the  joint  angle  6  and 
angular  velocity  9,  respectively.  The  affine  nonlinear  control 
system  is  realized  in  coordinates  x  =  (xi,  X2)  =  ( q ,  q)  by 


fix)  = 


x2 


G  .  .  .  k  b 


(54) 


g{x)  = 


i  L2 


(55) 


and  h  (x)  =  x±. 

The  Hamiltonian  H  for  the  pendulum  system  is  given  by 

H{x)  =  K  (x)  +  U  (x) 

=  ^mL2x\  +  ^  kx\  —  mGTcos(xi)(56) 

Lurthermore,  the  equipartition  of  energy  condition  (32)  is 
satisfied  trivially  for  the  1-DOL  system  with  ratio  £  =  b.  Ap¬ 
plying  Theorem  2.5,  the  controllability  function  Lc  is  given, 
exactly,  by 

Lc(x)  =  —  2bmGLcos(xi)  +  bkx\  (57) 
+b mL2  x 2  +2  bmG L 

Numerical  studies  were  conducted  via  application  of  our 
algorithms  to  the  pendulum  system  with  parameter  values 
given  in  Ligure  2.  The  observability  function  La  was  com¬ 
puted  using  an  algorithm  based  directly  on  Equation  (4).  The 
controllability  and  observability  functions  for  the  pendulum 
system  with  position  read-out  are  shown  in  Ligure  3. 

We  also  computed  an  approximation  to  the  controllability 
function  via  Equation  (15)  and  a  Monte-Carlo  approach.  We 
simulated  50,000  sample  paths  for  the  pendulum  system  with 
approximate  Gaussian  white  noise  injected  as  the  torque  in¬ 
put.  We  assumed  that  steady-state  was  reached  after  60  time 
units,  6  times  the  largest  time  constant  of  the  system. 

The  results  of  the  Monte-Carlo  experiments  are  presented 
in  Ligure  4.  In  this  case,  we  used  a  relatively  coarse 
grid,  which  is  roughly  the  highest  resolution  that  provides 
a  smooth  approximation.  By  generating  additional  sample 
paths,  we  can  increase  the  grid  resolution  while  maintaining 
a  smooth  approximation. 

Even  though  we  have  an  exact  expression  for  Lc,  we  use 
the  approximation  L'c  generated  by  the  Monte-Carlo  experi¬ 
ments  in  order  to  compute  a  balanced  realization  for  the  pen¬ 
dulum  system.  We  do  this  to  demonstrate  that  L'c  can  be  a 


Pendulum:  Exact  Controllability  Function 


Pendulum  (Position  Output):  Singular  Value  Function  1 


Figure  5:  The  singular  value  functions  for  the  pendulum  sys- 
Figure  3:  The  exact  controllability  function  (top)  and  observ-  tem  with  position  output.  Left:  a\  ( x )  nearly  constant  0.367; 

ability  function  (bottom)  for  the  pendulum  system  with  posi-  Right:  a>  ( x )  nearly  constant  0.284. 

tion  read-out. 


Stochastically  Excited  Pendulum:  Stationary  Density 


Pendulum:  Approximate  Controllability  Function 


Figure  4:  The  stationary  density  (top)  and  derived  approxi¬ 
mate  controllability  function  (bottom)  for  the  pendulum  sys¬ 
tem. 


suitable  surrogate  for  Lc  in  the  balancing  algorithms.  The 
algorithm  presented  in  Section  3  was  applied  to  L'c  to  com¬ 
pute  a  Morse  coordinate  transformation.  Further  computa¬ 
tions  (see  [3])  yielded  a  balanced  representation  for  the  pen¬ 
dulum  system  with  position  read-out. 

The  computed  singular  value  functions  (see  [1,  2])  are 
shown  in  Figure  5.  Because  the  pendulum  system  is  nearly 
linear,  we  expect  the  singular  value  functions  to  be  nearly 
constant  at  the  value  of  the  corresponding  Flankel  singular 
values  of  the  linearized  system.  This  is  reflected  in  the  com¬ 
putations.  The  singular  value  functions  are  nearly  constant 
at  grid  points  close  to  the  origin,  taking  values  close  to  0.367 
and  0.284.  This  closely  matches  the  Flankel  singular  values 
of  the  linearized  system  which  are  0.3671  and  0.2838.  One 
state  is  roughly  1 .3  times  as  important  to  the  input-to-output 
behavior  of  the  system. 

We  simulated  the  pendulum  system  in  the  original  and  bal¬ 
anced  coordinates  using  two  input  signals:  u  =  0  (natural 
response)  and  u(t)  =  0.5  sin  (t/n).  The  output  responses 
are  shown  in  Figure  6.  Theoretically,  the  output  responses  of 
the  original  and  balanced  systems  should  be  identical,  since 
they  are  merely  different  representations  of  the  same  physi¬ 
cal  system.  Flowever,  the  computations  introduce  numerical 
error. 

In  [3],  it  was  demonstrated  that  by  using  the  exact  con¬ 
trollability  function,  the  output  responses  of  the  original  and 
balanced  systems  are  virtually  identical.  Thus,  the  algo¬ 
rithms  for  computing  the  Morse,  input-normal,  and  balanc- 


Pendulum:  Output  (Position)  Response  -  Zero  Input 


Figure  6:  Output  response  for  the  pendulum  system  with  po¬ 
sition  read-out:  original  coordinates  (solid)  vs.  balanced  co¬ 
ordinates  (dashed).  Top:  zero  input,  approximate  Lc\  Bot¬ 
tom:  sinusoidal  input,  approximate  Lc. 


Figure  7:  Planar  double  pendulum  system  with  massless 
shafts,  linear  torsional  damping,  linear  torsional  stiffness, 
and  torque  input  applied  at  the  rotary  joints.  Parameters  val¬ 
ues  used  for  numerical  studies:  mi  =  m2  =  61  =  62  = 
fci  =k2  =  L  1=L2  =  1,G  =  10. 


The  double  pendulum  system  is  not  integrable  and  does 
not,  in  general,  satisfy  the  equipartition  of  energy  condition. 
However,  in  the  special  case  where  bi  =  b  =  b2  the  equipar¬ 
tition  of  energy  condition  is  satisfied  with  ratio  i  =  b.  Apply¬ 
ing  Theorem  2.5,  the  controllability  function  for  the  double 
pendulum  system  is  given,  exactly,  by 


ing  transformations  introduced  negligible  error.  On  the  other 
hand,  when  using  the  approximate  controllability  function 
generated  using  Monte-Carlo  data,  the  output  responses  of 
the  original  and  balanced  systems  deviate  somewhat.  Thus, 
a  better  approximation  may  be  desirable,  which  can  be 
achieved  by  generating  additional  Monte-Carlo  data. 

We  now  consider  a  double  pendulum  system  as  illustrated 
in  Figure  7.  As  with  the  pendulum  system,  the  system  incor¬ 
porates  linear  torsional  damping,  linear  torsional  stiffness, 
and  torque  inputs  at  the  rotary  joints.  We  assume  that  the 
shafts  are  massless  and  that  the  pendulum  moves  only  in  the 
plane.  We  measure  the  horizontal  position  of  the  end-effector 
as  the  system  output  (a  nonlinear  function  of  the  state  vari¬ 
ables). 

As  before,  we  obtain  a  state-space  realization  (/,  g,  h)  for 
the  pendulum  system  via  Lagrangian  mechanics.  The  affine 
nonlinear  control  system  is  realized  in  coordinates 

x  =  (x±,x2,  x3,  x4)  =  (qi,q2,qi,q2)  by 


/  0) 
9{x) 


q 

— M-1  (q)  (1 C(q,q)  +  N(q,q )) 

0 

M-1  (q)  _ 


(58) 

(59) 


where  expressions  for  the  matrices  M,  C,  and  N  can  be 
found  in  [3]  and  h  ( x )  =  L\  sin  (qi)  +  L2  sin  (q\  +  q2). 


Lc(x)  =  (60) 

b  (mi  +  m2)  L\  x\  +  bm2  L\  (x3  +  x4)2  + 

2  b  m2  Li  L2  cos  ( x2 )  X3  (X3  +  X4)  + 

bk\  x\  +  k2x\—2b  (mi  +  m2)  G  L\  cos  (x\)  — 

2bm2G  L2  cos  (x\  +£2)  + 

2bG  ((mi +m2)  Li +m2L2)  (61) 


A  numerical  study  was  conducted  via  application  of  our 
algorithms  to  the  double  pendulum  system  with  parame¬ 
ter  values  given  in  Figure  7.  Graphical  representations  of 
the  controllability  and  observability  functions  are  presented 
in  [3].  Figure  8  shows  a  2-dimensional  slice  of  each  sin¬ 
gular  value  function  for  the  double  pendulum  system.  At 
the  origin,  the  singular  value  functions  take  the  values,  re¬ 
spectively,  0.487,  0.444,  0.135,  and  0.050.  These  values  are 
reasonable  close  to  what  we  expect  from  the  Hankel  singular 
values  of  the  linearization,  which  are  0.5029, 0.4702,  0.0249, 
and  0.0106.  Two  states  of  the  balanced  realization  have  con¬ 
siderably  greater  input-to-output  importance  than  the  other 
two  states.  We  also  observe  that  numerical  errors  are  more 
prominent  for  the  singular  value  functions  of  small  magni¬ 
tude,  i.e.,  the  oscillations  that  they  display  are  likely  caused 
by  numerical  error  rather  than  being  an  accurate  reflection  of 
their  actual  behavior. 


Figure  8:  Singular  value  functions  for  double  pendulum  (2:3- 
X4  plane).  Top  left:  cr1(x);  Top  right:  <72(2:);  Bottom  left: 
03(2;);  Bottom  right:  0-4(2:). 

5  Concluding  Remarks 

Because  of  the  computational  complexity  of  our  algorithms, 
this  research  merely  represents  a  first  step  toward  making  the 
balancing  procedure  a  practical  model  reduction  tool.  It  is 
also  clear  that  results  with  broader  applicability  will  be  of 
benefit. 
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